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ABSTRACT 

Moment models of carrier transport, derived from the Boltzmann equation, have made 
possible the simulation of certain key effects through such realistic assumptions as energy 
dependent mobility functions. This type of global dependence permits the observation of 
velocity overshoot in the vicinity of device junctions, not discerned via classical drift-diffusion 
models, which are primarily local in nature. It has been found that a critical role is played 
in the hydrodynamic model by the heat conduction term. When ignored, the overshoot is 
inappropriately damped. When the standard choice of the Wiedemann- Franz law is made for 
the conductivity, spurious overshoot is observed. Agreement with Monte-Carlo simulation in 
this regime has required empirical modification of this law, as observed by IBM researchers, 
or nonstandard choices. In this paper, simulations of the hydrodynamic model in one and 
two dimensions, as well as simulations of a newly developed energy model, the RT model, will 
be presented. The RT model, intermediate between the hydrodynamic and drift-diffusion 
model, was developed at the University of Illinois to eliminate the parabolic energy band and 
Maxwellian distribution assumptions, and to reduce the spurious overshoot with physically 
consistent assumptions. The algorithms employed for both models are the essentially non- 
oscillatory shock capturing algorithms, developed at UCLA during the last decade. Some 
mathematical results will be presented, and contrasted with the highly developed state of 
the drift-diffusion model. 
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1 Introduction 


During the last decade, device modeling has attempted to incorporate general carrier heating, 
velocity overshoot, and various small device features into carrier simulation. The popular 
wisdom emerging from such concentrated study holds that global dependence of critical 
quantities, such as mobilities, on energy and/or temperature, is essential if such phenomena 
are to be modeled adequately. In this paper, we examine in detail the simulation of two 
such energy models, including the hydrodynamic model and the RT model. We describe the 
models, summarize some associated mathematical results, as well as the basic features of the 
numerical algorithm, and then present the results of extensive numerical simulations for two- 
dimensional MESFET devices, and for one-dimensional diodes. Both models represent one 
carrier flow. The hydrodynamic model contains hyperbolic modes related to the momentum 
equations, while the RT model does not possess such modes. In both cases, however, we 
employ a conservation law format, and numerical methods suitable for such systems. The 
ENO (essentially non-oscillatory) method employed makes use of adaptive stencils, and is 
particularly adept at shock capturing if the parameter regime crosses from supersonic to 
subsonic. Even if this does not occur, the convective terms are effectively discretized, via 
this procedure, in both models. The first use of such methods in device simulation was in 
[7], followed by the study [6], in which shocks were detected in micron devices at liquid Ni- 
trogen temperatures, and at room temperature in shorter devices, by independent numerical 
techniques. 

Our development of the RT model follows that of [5]. These researchers attempted 
to utilize a microscopic relaxation time approximation, which would allow for nonparabolic 
energy bands and non-Maxwellian distribution functions. The approach allows for parameter 
fitting of certain key quantities via Monte-Carlo simulation. 

One of the principal conclusions of the paper is the essential dependence of the hydrody- 
namic model upon the heat conduction term. Standard choices lead to numerically detected 
spurious overshoot at the drain junction of an n + — n — n + diode, while other choices signifi- 
cantly damp this overshoot. Monte-Carlo simulations show that substantial underestimation 
occurs when the heat conduction term is neglected. We refer the reader to [10], and to the 
simulation results of this paper for amplification. The RT model was developed, partly in 
response to the continuing debate concerning heat conduction processes in the hydrodynamic 
model. 

The status of mathematical results differs sharply between the hydrodynamic model, 
on the one hand, and the drift- diffusion model on the other. For the former, we have 
summarized two results, one by Gamba (cf. [8]) for an idealized model, in which the adi- 
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abatic relation is employed, and another by Gardner, Jerome, and Rose (cf. [9]) in which 
a Newton-Kantorovich theorem is developed for the n + — n — n + diode, yielding existence 
and convergence in a specialized subsonic regime. The drift-diffusion model, on the other 
hand, has been widely studied. Existence and approximation results have been carefully 
developed, although uniqueness is still not well understood for this model. Existence for 
the steady-state model is due in varying degrees of generality to many authors, including 
Mock ([17]), Seidman ([21]), and the first author ([13]). A convergence theory, based upon 
a calculus due to Krasnosel’skii, was presented in ([15]). Mathematical results have not yet 
been developed for the strongly nonlinear RT model. 

2 Hydrodynamic and Drift- Diffusion Models 

2.1 Mass, momentum and energy transport equations 

The equations as presented here are discussed in references [3], [20], and [4]. They are derived 
as the first three moments of the Boltzmann equation, with the latter written for electrons 
moving in an electric field as 

% + u-V x f--FVJ = C. (1) 

at m 

Here, / = f(x,u,t) is the numerical distribution function of a carrier species, x is the 
position vector, u is the species’ group velocity vector, F = F(x,t ) is the electric field, e 
is the electron charge modulus, m is the effective electron mass, and C is the time rate of 
change of / due to collisions. In the Boltzmann equation above, it has been assumed that the 
traditional Lorentz force field does not have a component induced by an external magnetic 
field. The moment equations, which will be derived subsequently, are expressed in terms of 
certain dependent variables, where n is the electron concentration, v is the average velocity, 
p is the momentum density, P is the symmetric pressure tensor, q is the the heat flux, ej is 
the internal energy, and C n , C p , and Cw represent moments of C, taken with respect to the 
functions 

h 0 (u) = 1, 
hi(u) = mu, 
h 2 (u) = 

The equations are given by: 

dn _ , . 

— + V- (nv) = C n , 


(2) 



d ,mn . , _ . ,mn . , . 

4 di^~ I V * +mne/ ) + V ‘ ( u I v I + mne /l) = 

— env • F — V- (yP) — V- q + Cjy. (4) 

The first Maxwell equation for the electric potential must be adjoined; each species con- 
tributes a corresponding moment subsystem, with appropriately signed charge. We begin the 
derivation with the definitions and assumptions. The concentration is given by n := f f du ; 
the average velocity by v := £ f uf du; the momentum by p := mnv ; the random velocity 
by c := u — v\ the pressure tensor by PiJ ■■= m / C{Cjf du; and the internal energy density by 
e/ := ^ / | c | 2 / du. This function represents energy/unit mass/unit concentration. The 
heat flux q is given by ^ := y /c, | e | 2 / du. Finally, for reference in subsequent subsec- 
tions, the electron current density is given by J := — env , and the energy flux is given by 

5 := /«{y | u | 2 }/ du. The assumptions on / are now stated. The function / is assumed 
to decrease sufficiently rapidly at infinity: 

lim hi(u)f(u) = 0, i = 0,1,2. 

|u|— KX) 

The derivation of (2), (3), (4) proceeds by multiplying the Boltzmann equation (1) by 
h 0 , hi, and h 2 , respectively, and integrating over group velocity space. With the application 
of certain standard identities ([16]), the mass/momentum/energy system is obtained. In 
addition to these transport equations, we have Poisson’s equation for the electric field, where 
rid := doping and t := dielectric: 

F = -W, (5) 

V-(eV^) = “ n d- (6) 

Here, we have used the convention that there are different species, each of concentration n,- 
and charge e,-. The entire system consists of equations (2), (3), (4), repeated according to 
species, and (5), (6). 

2.2 Moment closure and relaxation relations 

The system derived in the preceding subsection has fifteen dependent variables in the case 
of one species, determined by (j), n, v , P, e/, and q. By moment closure is meant the 
selection of compatible relations among these variables, so that the number of equations is 
equal in number to the remaining primitive variables selected. The relations to follow are 
characterized by the isotropic/parabolic energy band assumption. We begin by introducing 
a new tensor variable T, the carrier temperature, defined by 

Pij = nkTij 
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where k is Boltzmann’s constant, and a scalar variable W, the total carrier energy. A 
program of reduction to a set of basic variables, n, v, W, and <f>, or a set equivalent to these, 
can be implemented' by the following assumptions: 

• The pressure tensor is isotropic, with diagonal entries P a and off-diagonal entries zero, 
for a suitable scalar function, P s . P 3 is related to ej via mne/ = | P a . 

• It follows from the previous assumption that temperature may be represented by a 
scalar quantity T, and that the internal energy is represented in terms of T by 

met = -kT. 

2 

• The total energy density (per unit concentration) w is given by combining internal 
energy and parabolic energy bands with m assumed constant: 

w = me i -f -m | v | 2 , 
z 

and the total energy (per unit volume) W is the product, W = nw. 

• The heat flux is obtained by a differential expression involving the temperature: 

q = — kVT. 

Here, k is the thermal conductivity governed by the Wiedemann-Franz law (cf. [2]), 
described by 

K = ( _ + r)n _ T( _ ) . (7) 
The standard choice for r is r = — 1, but this has some associated difficulties. This 
will be amplified later in the paper. Here we simply remark that the term raised to 
the exponent r in (7) is proportional to the mobility, which in turn is proportional to 
the momentum relaxation time. 


In the case of N species, the closure relations determine ( d + 2 )N + 1 variables in d 
spatial dimensions. It is possible to rewrite the system (2, 3, 4) with the closure assumptions 


incorporated. We have the following. 

at + v ' 

dv 

-£ + v(V-p) + (p-V)v 
dw 

— + V- (u W) 


= c n . 

(8) 

= —enF — 'V(nkT) + C p , 

(9) 

= -env- F — V- ( vnkT ) 


+V- (kVT) + C w . 

(10) 
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The final step deals with the replacement of the collision moments. Motivated by the 
approach of [18], [1], [20], and [11], we define the recombination rate R and the momentum 
and energy relaxation times, r p and t w , respectively, in terms of averaged collision moments 
as follows. 



2.3 Drift-diffusion model 

The drift- diffusion model may be obtained by taking zeroth order moments of the BTE 
and adjoining the Poisson equation. Thus, one obtains the system for N carriers with 


concentrations n; and charge e,-, i = 1, • • • , N: 

-§f+ v 'J< =' -Ri, (13) 

F = -V*, (14) 

V-(eV« = -J2eim-n d . (15) 

There still remains the issue of determining the constitutive current relations. Classical 
drift-diffusion theory gives, for N = 2, n\ = n, and n 2 = p, 

J n = —ep n nV(f) + efitjVn, (16) 

J p = -ep,ppS/<f) - eDpVp. (17) 
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The introduction of exponential relations for n and p is also common, as is the use of the 
Einstein relations linking the mobilities, p n , p p , and the diffusion coefficients D n ,D p . These 


relations are specified by 



- 

D n = ( kT/e)p n , 

(18) 


D p = ( kT/e)p p . 

(19) 


It is also possible to derive the constitutive relations (16), (17), from the first order moment 
relations under the assumption that the momentum relaxation times tend to zero. The 
details are given in [20]. In fact, the constitutive relations include a heat flux term as well, 
which is suppressed at constant temperature. If it is not suppressed, one has an energy drift- 
diffusion model. In this derivation, one uses the definition of mobility in terms of relaxation 
time. 


3 RT Models 

In this section, we shall employ a microscopic assumption upon the momentum relaxation 
time, viz. , we shall assume that the collision term C in (1) is of the form, 

C = (20) 

T p 

where fx is the odd part of /. Note that this contrasts with the macroscopic assumption on 
t p , employed in the hydrodynamic model as described in Section 2. There, the representation 
defining r p was a post averaged expression. Here, the expression is employed in the averaging. 
In this case, we may obtain an expression for the energy flux S: 

S = -[nfi E F + V{nD E )], (21) 

where p E and D E are tensor expressions for mobility and diffusion, defined in terms of 
moments, and E represents average energy per unit concentration. The details are furnished 
in [5] . It is also shown there that the current density has the usual drift-diffusion form, with 
tensor expressions for mobility and diffusion. The RT model makes the following microscopic 
assumptions, with distinction between £ and its average, E. 

1. The even part of / is isotropic, and a function of £ alone, and the relaxation time is 
an inverse power function of £: 

fo = Me), 

r p = t p (£) = C£~ r . 
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2. The microscopic kinetic energy is a quadratic function of €, and the mass is not assumed 


constant: 


G(S) :=^M 2 =£ + «£ 2 , 


where a is an appropriate fitting parameter. 


3. The temperature is a modified variable in terms of which the following constitutive 
relation holds for E: 

E = 1 + h -akT). (23) 


Equation (23) allows for nonparabolic energy bands as well as non- Maxwellian distributions. 
Altogether, the model may be written in terms of the Poisson equation, (6), in conjunction 
with the system, 

V.J = 0, (24) 

V-S = U). (25) 

Here, J and S have been described previously, the latter in (21). In the expressions for J 
and S, the assumptions made for the model lead to scalar representations for the mobility 
and diffusion coefficients. For example, the choice made in [5] leads to 

H = /i 0 To/T, (26) 

/** = \(l-\kT)»kT. (27) 

The diffusion coefficients are defined by Einstein’s relations. The collision term in (25) is 

a specified quadratic function of T. One significant advantage of the microscopic relax- 
ation time (RT) assumption is that certain key parameters may be fitted via Monte-Carlo 
preprocessing, ensuring reliability of their values. 


4 Mathematical Results for the Hydrodynamic Model 

In this section, we shall describe some recent mathematical results, obtained in one spatial 
dimension. In the first subsection, we shall present existence and boundary layer results for 
a simplified version of the steady state hydrodynamic model. This will be followed by a 
convergence analysis for Newton’s method in the subsonic case. 

4.1 Existence and boundary layer theory 

We first write down the one dimensional evolution system in the case of a single carrier, in 
the absence of recombination. 


dn d(nv) 
dt dx 


= 0 , 
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(29) 


dp d(pv + knT) 
dt dx 

dW d(vW + vknT) 
8t + dx 

dF_ 

£ dx 


—enF — 

UX Ty; 

—en — rid • 


(30) 

(31) 


The corresponding steady state system is obtained by setting the time derivatives equal to 
zero. 

If we rewrite the second steady state equation by use of the pressure, P, we obtain the 
equation, 


d(pv -h P) _ 


dx 


= -enF 


(32) 


The approach of [8] is to eliminate the energy equation (30) from the system, and replace 
its role by a relation in the spirit of gas dynamics, i. e. by the constitutive relation, 


P(n) = ErP, 7 > 1. (33) 

When units are selected in which e = l,e=l,m = l, and K = 1, we obtain the system, in 
which n and <f> are the only dependent variables, 


j 

(TXn))* := (^ + n 7 )* 

<t>XX 


nv = constant, 

(34) 

j 

nfa S(fa, n), , 

(35) 

t p 

n — nj. 

(36) 


One can nominally specify boundary conditions on n and <j) at the endpoints of the device, 
taken here as the interval, [0, 1]. If the doping is such that the built-in potential is the same 
at both ends of the device, then we may take, 


*(0) = 0, *(1) = .ft, 


(37) 


for </> and 

n(0) = n 0 , n(l) = ni, (38) 

for n, where fa must satisfy the following consistency relation with respect to j , no, and n\\ 

iX (39) 


Here, 


fe = /(».,>)-/(»0,i)+i( ^,^ . 


f(nj) = ^2+7 n ' 1 2 - 


(40) 
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It is shown in [8] that a weak solution exists for the system (34), (35), (36), satisfying 
the boundary conditions exactly, or, in lieu of this, satisfying precise limiting relationships. 
The result for <j) I s classical, because the equation is elliptic. The result for n is provisional, 
and is detailed now. If j > 0, there is a weak solution n such that the relation, 

n = G + a, ' (41) 

holds, where G is Holder (|) continuous, and a is monotone increasing. Although n need 
not be of bounded variation, it is an entropy solution, in the sense that the function of x , 

H(n(x )) = ( F(n(x )) - F(n min ))signurn(n(x) - n min ) + Cx, (42) 


is monotone increasing. Here, n mtn is a minimum (location) for F , and C — sup S over 
relevant arguments. The following precise statement is available for subsonic boundary 
conditions. If no,ni > n m ,„, then the following holds. 


• Either 


n(l) = m, 


• or 


lim n(x) 

x-*l~ 

exists, and is a supersonic value, i. e. , is less than rc m ,- n , and even less than the conjugate 
value of n \ . 


• Similarly, either 


n(0) = n 0 , 


•. or 

lim n(x ) 

*-o+ v ' 

exists and is not less than n m ; n . 

In the second instance of both cases above, boundary layers occur, the one on the right 
involving transition through the supersonic regime. 

4.2 Newton convergence theory 

In this subsection, we order the basic variables as v, n, T, and <f>, because of symmetry 
considerations. Dirichlet boundary conditions are imposed, on n, T, <f>, with n(x m { n ) = 
n(x max ). In one dimension, the steady state equations are obtained from (28), (29), (30), and 
(31) by setting the time derivatives equal to zero. Dirichlet boundary conditions are imposed 
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in this subsection on n, T, and <f>, with n(x m i n ) = n{x max ). Since, by the conservation of 
mass equation, nv = j, it follows that the boundary conditions for v are periodic. The map 
defined by bringing all terms to the left hand side of the steady state system is called 4>. The 
linearized equations thus assume the form, where the boundary conditions are homogeneous 
Dirichlet conditions for 8n , 8T, and <5^>, and the boundary conditions on 6v are prescribed 
to be periodic, 


0 

0 

n ax ' ax ' 

+ 

A B 
c D 

d 

dx 

8v 

6n 

8T 

+ 

E F 
G H 


8v 

8n 

8T 

_ e 4f M 
t dx* J 




. W . 




. H . 


(43) 


The (spatially) dependent eigenvalues of the symmetric matrix A are calculated to be 


X=\(n + ~)± W(n + — Y 

2 V mnJ 2V \ mnj \rn J 


(44) 


Here, 



and the smaller eigenvalue is positive if n and T are strictly positive, and if 


(45) 


v 2 < — = c 2 . (46) 

m 

This type of point in function space is termed a subsonic point. This case was first considered 
in [9], where damped Newton/standard finite difference methods were presented. When 
Newton’s method is employed in this way, it is essential to determine conditions under 
which the linear increments are appropriately bounded. This is equivalent to uniform bounds 
for the operator derivative inverse maps, and represents one of the three properties for an 
(exact) operator Newton method to yield existence of a root and jR-quadratic convergence. 
The other two are sufficient regularity, and a sufficiently small starting residual, as measured 
in the range space norm. Explicit representations of B, C, D and of F, G, H are given in 
[14]. Moreover, if the system map $, subject to appropriate Dirichlet boundary conditions 
on n, T, and <f>, and periodic boundary conditions on v, accordingly has the domain D$ C 
X = Eli W 1,0 ° x fli W 2 '°°, and range in Y = IIi L °° , then (46) will hold for every element 
in a closed ball B ro C X , centered at a subsonic element uq € X, such that n and T are 
strictly positive, if r$ is sufficiently small. It is appropriate to assume at the outset, then, 
that D<t> C B r o, so that every function point in D<j> satisfies (46); we may also assume that 
n and T are uniformly bounded away from zero in this set. 

The Lipschitz property of the map <f>', where here we view the system (43) as the repre- 
sentation for §'(v,n,T,<f>)(z,u>) ~ /, with 


z = (Sv,Sn), u=(STM), / = (/i.A). (47) 
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is evident from the representation for 

The uniform inverse bounding proceeds as follows. As shown in [9], with the function 
spaces selected in this paper, the H 1 product norm of z can be estimated in terms of the 
L 2 norms of u>, u>' , and /i, under the conjunction of the hypothesis (46) and the L 2 x L 2 
coerciveness assumption, 

dA 

A = E + E* — is uniformly positive definite. (48) 

dx 

Here, E* is the matrix transpose of E and the latter is defined by 


E = 


dn 

dx 

dv | 1 

. dx ' Tp 



( 49 ) 


A final calculation, making use of the inner product of [0,u>] with (43), and the hypothesis, 


fin is positive, and sufficiently large, 


( 50 ) 


where h n is the nonzero entry of H, given by 


, dv ^r„-(W-W„K 

/in - — H o 


dx 


( 51 ) 


with t' w = shows that the H 1 product norm of [z, u] is estimated in terms of the L 2 norm 
of /. This series of calculations controls the L°° norm of [z,u>]; the L°° norm of [z 1 , to"] is 
now estimated by direct use of the system (43), making use of the fact that [v, n, T, <j> ] € B ro . 
We have now outlined the proof of the following. 


Theorem 4.1 Let the function spaces X and Y be selected as above, let the steady state 
system map $ be given with Dirichlet boundary conditions on n, T, (j>, and periodic boundary 
conditions on v, such that D$ C B ro , where every point in B ro is a subsonic point, with uni- 
form positivity bounds. If (48) and (51) hold, and x 0 is such that $(xo) is sufficiently small, 
then an R-quadratically convergent Newton sequence {i„} may be defined in the standard 
way, with limit x, satisfying $(x) = 0. 
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5 Discrete Schemes Based on Adaptive Stencils: ENO 


In this section, we shall briefly describe the ENO schemes as developed in [24] and [25]. 
Consider a system of hyperbolic conservation laws of the form 

d 

«* + £•/*■(“)*< = 9 (u,x,t), (52) 

i 


where 

u= x = ( ), 

and the hyperbolicity condition, 

d dfi 

y, is diagonalizable, with real eigenvalues, 
j ou 


holds for any real £ = (£j, • • • ,£*). An initial condition is adjoined to (52). 

For systems of conservation laws, local field by field decomposition is used, to resolve 
waves in different characteristic directions. Analytical expressions are employed for the 
eigenvalues and eigenvectors of an averaged Jacobian matrix. Typically, the Roe average 
[19] is employed. One feature of the ENO schemes in [24] and [25], which is distinct from 
the original ENO schemes of Harten et al [12], is that multidimensional regions are treated 
dimension by dimension: when computing fi(u) Xi in any particular direction, variables in all 
other directions are kept constant, and the Jacobians are treated in this direction. This, in 
essence, reduces the determination of the scheme to the case of a single conservation law in 
one spatial dimension. Thus, to describe the schemes, consider the scalar one dimensional 
problem, and a conservative approximation of the spatial operator given by 


Here, the numerical flux f is assumed consistent: 


f j+ L =/(«,■-/, •••,«,•+*); /(«,•••,«) = /(«)• (54) 

The conservative scheme (53), which characterizes the / divided difference as an approx- 
imation to f(u) x , suggests that / can be identified with an appropriate function h satisfying 

/(«(*)) = f + <J Mf) <*£, (55) 

If H is any primitive of h, then h can be computed from H' . H itself can be constructed 
by Newton’s divided difference method, beginning with differences of order one, since the 
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constant term is arbitrary. The necessary divided differences of H , of a given order, are 
expressed as constant multiples of those of / of order one lower. After the polynomial Q of 
degree r + 1 has been constructed, set 

= }• < 56 ) 
to obtain an rth order method. The construction is based on an adaptive stencil in the 
following sense: 

• One begins with an appropriate starting point to the left or right of the current “cell” 
by means of upwinding as determined by the sign of the derivative of a selected flux. 

• As the order of the divided differences is increased, the divided differences themselves 
determine the stencil: the “smaller” divided difference is chosen from two possible 
choices at each stage, ensuring a smoothest fit. 

• Lax-Friedrichs building blocks or Roe building blocks can both be used. For the lat- 
ter, in cells with sonic points, a local Lax-Friedrichs building block is used to avoid 
expansion shocks. 

Steady states are reached by explicit time stepping of arbitrary order; nonstandard high 
order Runge-Kutta methods exist [24] which preserve nonlinear stability of the first order 
Euler forward version under suitable CFL time step restrictions. The computer program 
is fully vectorized for computations on Cray supercomputers. For details of the efficient 
implementation, see [23]. 

6 Conservation Law Format for Hydrodynamic and 
RT Models 

In this section, we shall specify the conservation law format for the two dimensional hydro- 
dynamic model, and for the one dimensional RT model. 

6.1 Hydrodynamic model conservation format 

Define the vector of dependent variables as 

u = {n,<r,T,W), (57) 

where p = (<r, r). The system (8), (9), (10) can be written in the concise form, in two 
dimensions, as 

u t -I- fi(u) x + Mu) y = c(u ) + G(u, (j>) + (0, 0, 0, V- («VT)). (58) 
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The following identifications have been made in (58). 


ft \-t— — 5aW - eLLiI' i 

1 m’ 3 mn 2 mn ’ mn 3mn G 3 m 2 n 2 ’ 

2 m’ mn ’ 3 mn 2mn ’ 3mn 3 m 2 n 2 ’ 


(59) 

(60) 



IT - W 0 


£7(u) = (0, — enFi, — enF 2 , — enF • u). 


(61) 

(62) 


The eigenvalues and eigenvectors of f[ and f 2 are known (cf. [23]), and are readily 
incorporated into the field by field' decomposition required for the implementation of ENO. 


6.2 RT conservation format 

We shall present the conservation law form of the RT model. We begin with the vector form, 

ut + f(u) x = g{u) xx + h(u). (63) 

In equation (63), 


nE 

u = (en, ), 

TO 


f(u) = <!>'n (e/<(£), V E (E) + D(E)), 
g(u) = ( nD(E ), nD E (E)), 


h(u ) 


dE 


(0, enn(E)(<f>') 2 + -(n - n d )nD(E ) - n( | co //)). 


e ' ' ’ ' ’ dt 


(64) 

(65) 

( 66 ) 

(67) 


It can be shown that the left hand side defines a hyperbolic system, since the eigenvalues of 
f{u ) are real, for all positive n and T. 


7 Numerical Simulation Results 


We now present numerical simulation results for one carrier, two dimensional MESFET 
devices and one dimensional diodes. The third order ENO shock-capturing algorithm with 
Lax-Friedrichs building blocks, as described briefly in Section 5 and in more detail in [25], is 
applied to the hyperbolic part (the left hand side) of Equations (6.2) and (6.7). A nonlinearly 
stable third order Runge-Kutta time discretization [24] is used for the time evolution towards 
steady states. The forcing terms on the right hand side of (6.2) and (6.7) are treated in a 
time consistent way in the Runge-Kutta time stepping. The double derivative terms on the 
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right hand side of (6.2) and (6.7) are approximated by standard central differences owing to 
their dissipative nature. The Poisson equation (2.6) is solved by direct Gauss elimination for 
one spatial dimension and by Successive Over- Relaxation (SOR) or the Conjugate Gradient 
(CG) method for two spatial dimensions. Initial conditions are chosen as n = rid for the 
concentration, T = Tq for the temperature, and u = v = 0 (two spatial dimensions) or u = 0 
(one spatial dimension) for the velocities. A continuation method is used to reach the steady 
state: the voltage bias is taken initially as zero and is gradually increased to the required 
value, with the steady state solution of a lower biased case used as the initial condition for 
a higher one. 

7.1 Two dimensional MESFET 

We simulate, using the Hydrodynamic model (6.2)-(2.6), a two dimensional MESFET of the 
size 0.6 x 0.2pm 2 . The source and the drain each occupies 0.1pm at the upper left and the 
upper right, respectively, with the gate occupying 0.2pm at the upper middle (Figure 1, left). 
The doping is defined by rid = 3 x 10 17 cm -3 in [0, 0.1] x [0.15, 0.2] and in [0.5, 0.6] x [0.15, 0.2], 
and n d — 1 x 10 17 cm~ 3 elsewhere, with abrupt junctions (Figure 1, right). A uniform grid 
of 96 x 32 points is used. Notice that even if we may not have shocks in the solution, the 
initial condition n = rid is discontinuous, and the final steady state solution has a sharp 
transition around the junction. With the relatively coarse grid we use, the non-oscillatory 
shock capturing feature of the ENO algorithm is essential for the stability of the numerical 
procedure. 



Figure 1: Two dimensional MESFET. Left: the geometry; Right: the doping nj 

♦ 

We apply, at the source and drain, a voltage bias vbias = 2V. The gate is a Schottky 
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contact, with a negative voltage bias vgate = — 0.8F and a very low concentration value 
n = 3.9 x 10 5 cm -3 obtained from Equation (5.1-19) of [22]. The lattice temperature is 
taken as To = 300 °K. The numerical boundary conditions are summarized as follows (where 
$o = ^r-ln with fa = 0.138 x 10 -4 , e = 0.1602, and n,- = 1.4 x 10 lo cm~ 3 in our units): 

• At the source (0 < x < 0.1, y = 0.2): $ = $o for the potential; n = 3 x 10 17 cm -3 
for the concentration; T = 300° K for the temperature; u = 0 gm/ps for the horizontal 
velocity; and Neumann boundary condition for the vertical velocity v (i.e. || = 0 
where ft is the normal direction of the boundary). 

• At the drain (0.5 < x < 0.6, y = 0.2): $ = $o + vbias = $ 0 + 2 for the potential; 
n = 3 x 10 17 cm -3 for the concentration; T = 300 °K for the temperature; u = 0 gm/ps 
for the horizontal velocity; and Neumann boundary condition for the vertical velocity 
v. 

• At the gate (0.2 < x < 0.4, y = 0.2): $ = $ 0 + vgate = $ 0 — 0.8 for the potential; 
n = 3.9 x 10 5 cm -3 for the concentration; T = 300 ° K for the temperature; u = 0 fxm/ps 
for the horizontal velocity; and Neumann boundary condition for the vertical velocity 
v. 

• At all other parts of the boundary (0.1 < x < 0.2, y — 0.2; 0.4 < x < 0.5, y — 0.2; 
x — 0,0 < y < 0.2; x = 0.6,0 < y < 0.2; and 0 < x < 0.6, y = 0), all variables are 
equipped with Newmann boundary conditions. 

The boundary conditions chosen are based upon physical and numerical considerations. 
They may not be adequate mathematically, as is evident from some serious boundary layers 
observable in Figures 2 through 6. ENO methods, owing to their upwind nature, are robust 
to different boundary conditions (including over-specified boundary conditions) and do not 
exhibit numerical difficulties in the presence of such boundary layers, even with the extremely 
low concentration prescribed at the gate (around 10 -12 relative to the high doping). We point 
out, however, that boundary conditions affect the global solution significantly. We have also 
simulated the same problem with different boundary conditions, for example with Dirichlet 
boundary conditions everywhere for the temperature, or with Neumann boundary conditions 
for all variables except for the potential at the contacts. The numerical results (not shown 
in this paper) are noticeably different. This indicates the importance of studying adequate 
boundary conditions, from both a physical and a mathematical point of view. 

In Figures 2 through 6, we shpw pictures of the concentration n, temperature T, horizontal 
velocity u, vertical velocity v, and the potential $. Surfaces of the solution are shown at 
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the left, and cuts at y — 0.175, which cut through the middle of the high doping “blobs 
horizontally, are shown at the right. 



Figure 2: Two dimensional MESFET, concentration n. Left: surface of the solution; 
Right: cut at y = 0.175 





Figure 3: Two dimensional MESFET, temperature T. Left: surface of the solution; 
Right: cut at y = 0.175 




Figure 6: Two dimensional MESFET, potential Left: surface of the solution; Right: 
cut at y — 0.175 

Notice that there is a boundary layer for the concentration n at the dram but not at the 
source. Also notice the rapid drop of n at the depletion region near the gate. The temperature 
achieves its maximum around the left corner of the drain. The leakage current at the gate 
appears negligible from the normal velocity component, while the horizontal component 
shows evidence of strong carrier movement toward the source beneath the left gate area, and 
strong movement toward the drain immediately to the left of the drain junction. 

We have also simulated the same MESFET with a higher doping ratio: 3 x 10 17 cm 3 in 
the high doping region versus 1 x 10 15 cm~ 3 in the low doping region. We observe similar 
results (pictures not shown here). 

7.2 HD model for a one dimensional diode — spurious velocity 
overshoot 

A notorious phenomenon of HD models is that spurious velocity overshoot occurs at the 
drain junction of an n + -n-n + diode. It is intrinsic to the model and is not a numerical 
artifact, as is verified by our grid refinement study and by using different numerical algo- 
rithms. This phenomenon is closely related to the physical assumption governing the heat 
conduction term. Gnudi, Odeh and Rudan [10] observed that the spurious overshoot can be 
greatly reduced by an empirical modification of the Wiedemann-Franz law for the thermal 
conductivity. 

In this subsection we perform an extensive numerical study of the dependency of the 
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spurious velocity overshoot upon the heat conduction term. The n + -n-n + diode we simulate 
has a length 0.6^m, with a doping defined by rid = 3 x 10 17 cra -3 in [0, 0.1] and in [0.5, 0.6], and 
n d = 1 x 10 15 cm -3 in [0.15,0.45], with smooth junctions (Figure 7). The lattice temperature 
is taken as To = 296. 21°/^. We apply a voltage vbias = 1 .5 V, as is the case in [10]. 



Figure 7: The doping rid for the one dimensional n + -n-n + diode 
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Figure 8: HD for one dimensional n + -n-n + diode. Velocity u (upper), temperature T 
(middle left), concentration n (middle right), total energy W (lower left) and electric field 
— (lower right). Solid line: r=-l; dashed line: r=-2; circles: r=-2 in the coefficient of (2.7); 
pluses: r=-2 in (2.7) 







The standard HD model uses r = — 1 in the Wiedemann-Franz law (2.7) and the relax- 
ation times (2.11). The numerical solution of this is shown as solid lines in Figure 8. We 
can clearly observe the spurious velocity overshoot at the right junction, but otherwise the 
solution is basically correct comparing with direct Monto-Carlo simulations (not shown). 
When r is taken as —2 in (2.7) and (2.11), the solution is completely wrong (dashed line in 
Figures 8). However, when one takes r = — 2 only in the coefficient of k in (2.7) but leaves 
r = — 1 in the power of k in (2.7) and in (2.11), i.e., when one uses 

_ 1 T ,_ 1 . . 

2 e ^ To ’ (68) 

in the place of (2.7) and leaves r — — 1 in (2.11) unchanged, as was done in [10], one obtains a 
greatly reduced spurious overshoot (the circles in Figures 8). Finally, the result with r = —2 
in k in (2.7) but with r = — 1 in (2.11) unchanged, is shown by pluses in Figure 8. We can 
see that the spurious overshoot also disappears. 

7.3 RT model for a one dimensional diode 

We present numerical simulation results for the RT model, described in Section 3, for the 
same one dimensional diode used in Subsection 7.2. Although the RT model is a parabolic 
system with two equations, the existence of sharp transition regions near the junctions 
justifies the usage of ENO shock capturing algorithms for the hyperbolic part. 

In Figure 9, we show the results of velocity u, temperature T, concentration n, total 
energy E, and electric field — of the RT simulation, in circles, in a background of standard 
HD results ( r = —1) in solid lines, and of HD results with r = — 2 in the coefficient of k in 
(2.7) but with r = — 1 in the power (i.e., (2.7) is replaced by (7.1)), and r = — 1 in (2.11), in 
dashed lines. We can see that the RT model greatly reduces the spurious velocity overshoot 
and is comparable with the result of the empirically modified HD result in dashed lines. 

Extensive numerical tests about the RT model, as well as comparisons between the RT 
and HD models, constitute ongoing research, jointly with U. Ravaioli, E. Kan and D. Chen 
at the University of Illinois. 

Acknowledgements: We would like to thank Edwin Kan, Umberto Ravaioli and Stan- 
ley Osher for helpful discussions. Computations were performed on the Cray YMP at the 
Pittsburgh Supercomputing Center and on the Cray YMP at the University of Illinois. 


22 


Figure 9: One dimensional n + -n-n + diode. Velocity u (upper), temperature T (middle 
left), concentration n (middle right), total energy E (lower left), and electric field — (lower 
right). Solid line: standard HD with r=-l; circles: RT; dashed lines: HD with r=-2 in the 
coefficient of (2.7) 
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